Population genetic diversity and structure of Rhinogobius candidianus (Gobiidae) in Taiwan: Translocation and release

Abstract Rhinogobius candidianus is a freshwater goby distributed in north, northwest, west, and south Taiwan, but this species has been introduced to east Taiwan and became dominant. To investigate its native population genetic diversity and structure and evaluate the sources and diversity of translocated populations, the mitochondrial DNA control region and cytochrome b gene (1981 bp) from 220 specimens were analyzed. These results indicated that (1) the east populations originated from two sources in west Taiwan; (2) translocated populations exist in east Taiwan and south Taiwan; (3) many populations have likely been moved secondarily by human intervention; (4) the effective size of the populations had declined greatly; (5) within the native populations, the ancestral populations colonized Taiwan during the land bridge phase in the Pleistocene through north Taiwan; (6) the landform changes in Taiwan shaped the population structure; and (7) the landforms of the coastline during glaciation also shaped the native range. The low‐level genetic diversity, high population differentiation, and population decline greatly suggest the need for resource management and conservation interventions. Four clades (α–δ) should be managed as four distinct evolutionarily significant units, while the translocated populations should be managed as separate management units. Moreover, the translocated populations in east Taiwan should be evaluated and monitored carefully.


| INTRODUC TI ON
Rhinogobius is a genus of freshwater goby native to East Asia. Most species are primarily freshwater fishes, and very few species are amphidromous. In the Catalog of Fishes (Fricke et al., 2019), there are currently 83 valid species listed within this genus. In FishBase (Froese & Pauly, 2022), there are 108 scientific names, but only 65 are valid species. Currently, new Rhinogobius species are being identified (Takahashi & Okazaki, 2017;Xia et al., 2018), and these new species descriptions may suggest the list is not yet complete.
Of the 10 described Rhinogobius species on Taiwan island, five species have been analyzed for their phylogeographic relationships; R. maculafasciatus (Cheng et al., 2005a), R. rubromaculatus (Cheng et al., 2005b), R. giurinus , R. delicatus (Ju et al., 2021), and R. gigas (Liao et al., 2021). Recently, our study found that the Rhinogobius species have been exported for the worldwide aquarium trade from Taiwan (our observations; Yang et al., 2020). Besides R. giurinus, an amphidromous species, the distribution area of R.
candidianus is wider than that of other Taiwan species. Chen and Shao (1996) proposed that R. candidianus was distributed in north, northwest, west, and south (north of Tzengwen River) Taiwan, but Leander et al. (2014) found that this species was also distributed in east Taiwan, and in and south of the Tzengwen River. Liao et al. (2021) proposed that R. candidianus has been introduced to east Taiwan and became dominant. Thus, our study analyzed mitochondrial DNA in order to characterize the genetic diversity of R. candidianus and determine the geographic origin of translocated stocks as well as any human-mediated gene flow between translocated stocks and native populations.
Previous studies (Chang et al., 2016;Chiang et al., 2010Chiang et al., , 2013Chiu et al., 2017;Han et al., 2019) have suggested that Taiwan island provides an excellent opportunity for examining phylogeographic patterns. Taiwan island is located off the southeastern coast of mainland China and the southern East Asian islands. Geological evidence indicates that land bridges connected this island to the Asian continent and the Japanese islands during glaciations (Fairbanks, 1989;Gascoyne et al., 1979;Kimura, 2000;Ota, 1998). Past phylogeographic studies have highlighted the influence of natural geological barriers, and general heterogeneous topography of Taiwan, in shaping the regional biogeography of freshwater fishes (Chiang et al., 2010(Chiang et al., , 2013Chiu et al., 2017;Han et al., 2019;Ju et al., 2018;Lin et al., 2016). Many ichthyofaunal and phylogeographic studies have suggested that geological barriers such as the Central Range, Miaoli Plateau, and Kaoping foreland basins shaped the structures and distribution patterns (Han et al., 2019;Ju et al., 2018;Lin et al., 2016) (Figure 1). Thus, the distribution patterns of freshwater fishes in Taiwan vary. For example, Opsariichthys pachycephalus, Candidia barbatus, Cobitis sinensis, and Acrossocheilus paradoxus are F I G U R E 1 Sampling localities of the Rhinogobius candidianus in Taiwan. Black dots represent approximate sampling sites. Sample size is given in Table 1. The black thicker dashed lines represent the border of each zoogeographical district, which were defined by Tzeng (1986). The gray thinner dotted lines represent the coastline during the glaciations following Ota (1998). The organism photograph was provided by Mingtai Zhou. widely distributed in north, west, and south Taiwan; Sinogastromyzon puliensis is distributed south of the Miaoli Plateau; and O. evolans, Squalidus argentatus, Sinibrama macrops, and Hemibarbus labeo are only distributed in the Tamsui River, north of the Taoyuan Plateau. Lin et al. (2016) proposed that the different distribution patterns of freshwater fishes in Taiwan resulted from different colonization times and colonization routes contributed by the geological history.
Thus, our study attempted to examine the intraspecific genetic diversity of native R. candidianus to understand its native distribution pattern, structure, diversity, and colonization history. Accordingly, the present major aims were to determine (1) the native, translocated, and released populations of R. candidianus, (2) the demography and genetic diversity of native R. candidianus, and (3) the genetic structure and colonization history of native R. candidianus in Taiwan. To achieve the above aims, the genetic variability was analyzed using mitochondrial DNA (mtDNA) cytochrome b gene (cyt b), and control region (CR or d-loop). The regions of the mtDNA are often analyzed in studies of animal phylogeography (e.g., Chen et al., 2007;Hsu et al., 2014;Yang et al., 2016). Moreover, these two loci have previously been used in the analysis of Taiwanese freshwater fishes and freshwater gobies (Chiang et al., 2010(Chiang et al., , 2017Ju et al., 2018Ju et al., , 2021Liao et al., 2021).

| MATERIAL S AND ME THODS
The care and use of experimental animals complied with Taiwan's animal welfare laws, the Taiwan Council guidelines and policies as approved by MOST 106-2611-M-291-007. Taiwan's government does not require any license to work with dead animals from fisheries. All methods and the handling of fish were performed and approved by the National Museum of Marine Biology and Aquarium.

| Population sampling and molecular methods
A total of 220 specimens of R. candidianus were collected from 21 localities across Taiwan (Table 1, Figure 1). All specimens are housed in the laboratory of Yu-Min Ju, National Museum of Marine Biology and Aquarium. These specimens were classified into five groups, north, northwest, west, south, and east, based on previous ichthyofaunal and phylogeographic studies (e.g., Han et al., 2019;Tzeng, 1986;Wang et al., 2004)  and (21) Beinan, BN. Fish were collected from field sites with seines and fatally anesthetized with MS-222 (Sigma). The samples were fixed and stored in 100% ethanol. Genomic DNA was extracted from muscle tissue using a genomic DNA purification kit (Gentra Systems, Valencia, CA). The entire cyt b gene and CR region were amplified by polymerase chain reaction (PCR) using primers from Xiao et al. (2001) and Zhou et al. (2017). Each 50 μl PCR reaction mixture contained 5 ng of template DNA, 5 μl of 10x reaction buffer, 4 μl of dNTP mix (10 mM), 5 pmol of each primer, and 2 U of Taq polymerase (TaKaRa, Taq polymerase). The PCR was programmed on an MJ Thermal Cycler as 1 cycle of denaturation at 94°C for 4 min, 40 cycles of denaturation at 94°C for 30 s, annealing at 51°C-55°C for 50 s-1 min, and extension at 72°C for 1 min 30 s, followed by a 72°C extension for 10 min and 4°C for storage. The purified PCR products were sequenced using an ABI 377 automated sequencer (Applied Biosystems, Foster City, CA, USA). The chromatograms were assessed using the software CHROMAS (Technelysium), and the sequences were manually edited using BioEdit 6.0.7 (Hall, 1999).

| Sequence alignment and phylogenetic inferences
The nucleotide sequences were aligned with Clustal X 1.81 (Thompson et al., 1997). The most appropriate nucleotide substitution model was TrN + I + G using the Bayesian information criterion (BIC) for cyt b, CR, and concatenated data sets in jmodelTest 2.0 (Darriba et al., 2012). The phylogenetic analysis was performed using a maximum likelihood (ML) estimation with MEGA-X (Kumar et al., 2018). Bootstrapping was performed with 1000 replications.
The median-joining algorithm (Bandelt et al., 1999) from Network 5.0 was used to reconstruct the haplotype networks. Additionally, the phylogenetic tree was also generated by the program BEAST 1.8.0  with 10 7 MCMC steps and the first 10% as burn-in. The strict clock model with a Bayesian skyline tree was used to construct Bayesian skyline plots and estimate the divergence times of the major lineages to the most recent common ancestor (T MRCA ) by running 10 6 generations. The molecular clock was calibrated using a divergence rate of 1.33% per million years for the concatenated data set (Burridge et al., 2008;Ju et al., 2021).
The output was visualized in Tracer v1.6  to determine that convergence and suitable effective sample size were achieved for all parameters. The burn-in and plots for each analysis were visualized using Tracer. The TreeAnnotator in the BEAST package was used to summarize the tree data, and the tree was viewed using FigTree v1.3 (Rambaut, 2014).

| Population diversity and structure
The levels of intrapopulation genetic diversity were estimated based on indices of haplotype diversity (h) (Nei & Tajima, 1983) and nucleotide diversity (θ π ) (Jukes & Cantor, 1969) in DnaSP 4.10.8 (Rozas et al., 2003). Two genetic differentiation indices, G ST and N ST , were used to examine the existence of a phylogeographic structure in DnaSP (Pons & Petit, 1996). Pairwise F ST and p-distance values implemented by DnaSP were used to examine the spatial partitioning of genetic variation among populations. To determine the possible diversification scenarios, a statistical dispersal-vicariance analysis (S-DIVA) was employed to determine statistical support for ancestral range reconstructions (Yu et al., 2010). The S-DIVA analysis was implemented in RASP v.3.1 (Yu et al., 2015). Range information was defined using the ichthyofaunal classification (Tzeng, 1986) (Table 1, Figure 1). The analysis was performed using the 'maxareas = 4-7' option.

| Population demography and history
The past demographic expansions were examined using three different approaches. First, Tajima's D (Tajima, 1989) and Fu's Fs (Fu, 1997) tests were used in DnaSP. Second, the demographic history was examined using mismatch distribution analyses implemented in DnaSP.
Thirdly, the effective population size changes over time were evaluated using the Bayesian Skyline Plots (BSP) computed with BEAST.
A Bayesian Skyline tree was selected and a strict clock model was used. We ran the analysis for 3 × 10 7 generations to ensure convergence of all parameters (ESSs > 200) and the first 10% of samples for each chain were discarded as burn-in. We used the same setting for two independent runs to check for convergence with 3 × 10 7 and sampling every 1 × 10 3 generations. Plots for each analysis were drawn using Tracer.
The Approximate Bayesian Computation (ABC) framework was used to determine the population history with the software DIYABC ver. 2.0 (Cornuet et al., 2014). The reference table was built with 1,000,000 simulated data sets per scenario using all statistics. We and 840 bp of the CR (66 variable sites and 55 phylogenetically informative sites)] from 220 specimens were analyzed. In the cyt b data set, a total of 55 haplotypes were obtained (refer to Table 1).
The ML phylogenetic tree revealed that these 55 haplotypes fall into three lineages (lineages I, II, and III) ( Figure 2a). The Bayesian tree revealed the same topology (data not shown). These three lineages were sympatric in some populations ( Figure 3a). Among the 55 haplotypes, only nine were shared by two or more populations, and six shared haplotypes (S01, S03, S06, S07, S08, and S09) were only distributed in adjacent populations ( Figure 3c). The S02, S04, and S05 were distributed in the west and east populations. Five populations (SS, FA, WS, LB, and BN) did not have private haplotypes (Table 1).
Four populations (FA, WS, JS, and LB) only included one haplotype.
In the CR data set, a total of 78 haplotypes were obtained (refer to Table 1). The ML phylogenetic tree revealed that these haplotypes fall into two major lineages (CR1 and CR2), but not all clusters were highly supported by bootstrapping test (Figure 2b Table 1). The population LB included only one haplotype ( Table 1).
Based on the results obtained from the cyt b and CR data sets (Table 1 (Figure 3b). To identify the relationships among these populations, our study investigated genetic diversity and population structure further using the concatenated data set. In the mtDNA haplotype network, all 91 mt haplotypes were assorted as three clades (Figure 4). This pattern was the same as the pattern of cyt b data set (Figures 2a and 3a).  1 and 4).
Within clade C, the populations HL, TA, and TJ were not close to population TW in geography, but the haplotype network grouped these three populations close to population TW (Figures 1 and 4).
Within clade A, the population ML was not close to populations FA, TC, and CK in geography, but the haplotype network grouped ML with specimens from populations FA, TC, and CK (Figures 1 and 4).
The pairwise p-distance values among populations also supported these results within each of the clades ( Table 2). Although disjunct distributions could be the result of extinctions, our study considered that these specimens in population ML within clade A, and these specimens in populations HL, TA, and TJ within clade C were secondarily moved through human intervention. Moreover, within clade A, the results of the pairwise p-distance revealed that the populations PZ and BC were close to the populations LB, HD, and BN (Table 2), although the haplotype network showed that they might originate from population CS (Figure 4). The population CS only included the cytochrome b phylogenetic lineage I, but the populations PZ and BC included lineages I and III (Figures 2a and 3a). Thus, our study also suggests that some specimens in populations PZ and BC originated

| Genetic diversity and structure
Our study assorted all sampling specimens as three groups, native, released, and translocated populations ( Table 3). Removing the released and translocated populations, a total of 151 native specimens were re-analyzed using concatenated data ( Table 3). The haplotype network based on mtDNA data set revealed these specimens fell into four allopatric clades, α-δ ( Figure 5). The clades α-δ were distributed in north, northwest, west, and southwest regions, respectively ( Figure 5). The clades β and γ were clustered together as  (Table 4).
These results showed that the population differentiations were high.
In all samples based on the concatenated data, Tajima's D and Fu's Fs tests, and mismatch distribution displayed non-signature of recent demographic expansion ( Figure 6). Moreover, the pattern of the mismatch distribution was multimodal, and these results supported the population differentiation (Yan et al., 2020;Yi et al., 2021). Thus, our study analyzed the demographic history in four clades based on native, native + released, and all samples (including native, released, and translocated) (Heller et al., 2013). The BSP, which simulated the fluctuations in population size over time, concluded population declined patterns in the four clades (clades α-δ, Figures 5 and 7). These results revealed that in the past, the clade α had higher effective population size, and the effective population sizes of other three clades were similar. However, today, the clade δ had higher effective population size relative to the other clades. In clade γ, the analyses of the BSP suggest that the effective population size in the all samples' data set, including the translocated and released populations, revealed higher population size than that in the native + released TA B L E 2 Matrix of p-distance (*10 −2 ; below) with concatenated clades A, B, and C in Figure 4 Clade A   Chen and Shao (1996)  Rhinogobius delicatus is a primary freshwater fish species endemic from Taiwan (Ju et al., 2021); and R. giurinus is an amphidromous species . In addition, the overall values of F ST in many primary freshwater species were similar to R. delicatus (e.g., giurinus was distributed widely on the island , and that of R. delicatus was only shared between adjacent populations (Ju et al., 2021). Thus, our study considered that R. candidianus is a primary freshwater fish and Taiwan endemic species. Chen and Shao (1996) discovered R. candidianus, and proposed this species was not distributed in east Taiwan, and in and south of the Tzengwen River (TW). However, Leander et al. (2014) found that this species was distributed in Taiwan widely and Liao et al. (2021) proposed that R. candidianus has been introduced to east Taiwan and became dominant. Actually, previous phylogeographic studies also found that other freshwater fish in west Taiwan were introduced to east Taiwan (Zacco pachycephalus see Wang et al., 1999; Varicorhinus barbatulus see Wang et al., 2004; Acrossocheilus paradoxus see Ju et al., 2018). Besides freshwater fishes, human introduction of a freshwater prawn from west to east Taiwan was also found (Macrobrachium asperulum see Liu et al., 2011). Ju et al. (2018) suggested that populations secondarily moved through human contact would affect the interpretation of the phylogeographic history.

| Native, translocated, and released populations
Thus, our study attempted to find the native distribution pattern to infer the phylogeographic history of R. candidianus. Chen and Shao (1996) found that R. candidianus was distributed north of the Tzengwen River (population TW) and did not find R. and BC were defined as released populations because these specimens were distributed in the native area but showed discordant phylogenies when compared to other geographically close populations (Table 3).
Among six released populations, population HL only included released specimens, but the other populations contained native and released specimens (Table 3) (Figures 2 and 4). It is possible that the population HL may have undergone a bottleneck, and then expanded, but its history needs to be tested further.

| Phylogeography of R. candidianus
Previous studies found that the mutation rate of mtDNA for the molecular clock calibration was variable, ranging from 0.54% (Ruber et al., 2007;Zardoya & Doadrio, 1999) (Ju et al., 2021). The results of the estimated T MRCA based on the concatenated data set suggested that the origin of R.
candidianus could be traced back to late Pleistocene (T MRCA = 0.94 mya), agreeing within previous observations (e.g., Chang et al., 2016;Ju et al., 2016Ju et al., , 2018. Thus, our study considered that this molecular clock is likely a proper estimate. Taiwanese orogeny (mountain building) uplifted the longitudinal Central Range to nearly 4000 m at 2 mya (Teng, 1990). The distribution patterns of freshwater fishes and phylogeographic studies indicate that the Central Range has acted as a barrier to dispersal between the western and eastern populations of species (e.g., Chang et al., 2016;Tzeng, 1986 1 and 5). This distribution pattern also suggested that R.
candidianus colonized Taiwan after 2 mya. Besides the Central Range, the landform in the east coast of Taiwan island also broke the dispersions of R. candidianus from northeast to east because the land along the east coast was not exposed even during glaciations ( Figure 1).
In total, 151 specimens from 16 native populations fell into four allopatric clades (Table 3, Figure 5). According to all results of the mtDNA haplotype network and ABC approaches ( Figures 5 and 8 , 8 and 9). Moreover, the southward colonization route was broken by the Kaoping foreland basins (Figures 1 and 5). The Kaoping foreland basins, reaching 200 m depth within 3 km of shoreline, formed in 2-3 mya and could be a barrier for the southern part of the Kaoping river (Boggs et al., 1979;Chiang et al., 2013;Hsu et al., 2014;Ju et al., 2018;Lin et al., 2016;Wang et al., 2004). Previous studies have documented that the Kaoping foreland basins shaped the intraspecific structure and the distribution area of many freshwater species (e.g., Chiang et al., 2010Chiang et al., , 2017Han et al., 2019;Hsu et al., 2014;Lin et al., 2016).
After southward dispersions, the Taoyuan Plateau formed and the populations in north Taiwan (clade α) diverged ( Figures 5 and 9).
The Taoyuan Plateau is located in northwest Taiwan (Figure 1) Chang et al. (2016) and Hsu et al. (2014) also found that the Taoyuan Plateau divided the populations of M.
brevirostris and S. libertina into different lineages. The present study found that the Taoyuan Plateau isolated the R. candidianus populations in the Tamsui River as clade α ( Figure 5).
The Formosa Bank is located in the south Taiwan Strait. Previous studies suggested that the Formosa Bank divided the glacial land bridge in the Taiwan Strait, but the effect of the Formosa Bank on the population dispersion within the island has not been described (e.g., Chang et al., 2016;Chiang et al., 2010;Lin et al., 2016;Oshima, 1923). However, Ju et al. (2018) proposed that during the maximum glacial period, the ridge uplifted from the Formosa Bank to the present coastline of Taiwan island. During maximum glaciation, the dispersions between the two sides of the bank through the exposed continental shelves of the island were fragmented.
Thus, the populations south of the bank were isolated and diverged as clade δ (Figures 5 and 9). Finally, many studies suggest that the Miaoli Plateau isolated the dispersion of freshwater fishes (Chang et al., 2016;Jean et al., 2014;Lin et al., 2016). Thus, when the Miaoli Plateau emerged, the populations were isolated and began to diverge as clades β and γ (Figures 5 and 9).

| Demography and conservation
The results of the Bayesian skyline plots showed that the population sizes of R. candidianus declined greatly (Figure 7). The neutral and mismatch distribution tests also supported the pattern of population decline ( Figure 6). The existence of low genetic diversity, a pattern of population decline, and high population differentiation support the needs for the development of management strategies. Four clades in the haplotype network based on concatenated data ( Figure 5) indicate that these four clades should be recognized as four evolutionarily significant units (ESUs) for conservation (Moritz, 1994;Ryder, 1986).
Besides in situ conservation, released and translocated populations are also methods for resource management. Our study found that the released and translocated populations, for example, HL, HD, SK, and BN, have many private haplotypes (Tables 1 and 3). The results of phylogenetic analyses also displayed that these populations were different than other native populations (Figures 2a and   4). In addition, our study found that three shared haplotypes (S03, T03, and T05) were only distributed in the east populations (translocated populations HD, SK, and BN; Figure 3c). These three haplotypes may have recently diverged from a common source in the east populations, and these three east populations likely originated from population HL (Figure 4). The population HL did not have shared haplotypes with any population (Figure 3c and d), and the phylogenetic tree and haplotype network branches among populations HL and other native populations were longer (Figures 2 and 4) (Palsbøll et al., 2007).
In addition, our study suggested that the translocated populations need to be carefully evaluated. In clade δ, the analyses of BSP displayed that the results in the two data sets, including and excluding the translocated populations (all and native + released), revealed the similar effective population sizes (Figure 7d). However, in clade γ, the BSP results displayed that the results in the all samples, including translocated populations, had higher effective population sizes than that in native + released samples, excluding translocated populations ( Figure 7c). Our study found that the nucleotide diversities of the translocated populations in the clade γ were higher than those in the clade δ (Table 3). Conclusively, if the nucleotide diversities in sources of the translocated populations were low, the work of ex situ conservation is not necessarily (Figure 7c and d). Besides, our study considered that these translocated populations in east Taiwan also need to be carefully evaluated because our study found that the body sizes are larger than those in west Taiwan

ACK N OWLED G M ENTS
We thank Ms. Shun-Ya Jian for help with the sample collection and Mr. Mingtai Zhou for providing the organism photograph. This research was supported by grants from the Ministry of Science and Technology of Taiwan (No. MOST 106-2611-M-291-007). We also thank the associate editor and anonymous referees for their helpful comments.

CO N FLI C T O F I NTE R E S T
The authors have no conflicts of interest to declare.

This article has earned Open Data, Open Materials and Preregistered
Research Design badges. Data, materials and the preregistered design and analysis plan are available at GenBank (https://www.ncbi. nlm. nih.gov/) under the accession numbers OK375877-OK375931, and OK375932-OK376009.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data set, which included 55 cytochrome b and 78 d-loop sequences, was submitted to GenBank (https://www.ncbi.nlm. nih. gov/) under the accession numbers OK375877-OK375931, and OK375932-OK376009.